CRS
Diego Hernangómez
2022-01-20
Un sistema de referencia de coordenadas (o CRS por sus siglas en inglés, Coordinate Reference System) permite relacionar datos espaciales con su localización en la superficie terrestre. Los CRS constituyen por tanto un aspecto fundamental en el análisis y representación de datos espaciales, ya que nos permiten identificar con exactitud la posición de los datos sobre el globo terráqueo.
Así mismo, cuando se trabaja con datos espaciales provenientes de distintas fuentes de información, es necesario comprobar que dichos datos se encuentran definidos en el mismo CRS:
Figure 1: Representación de mismos valores de coordenadas en distintos CRS
En la Fig. 1, ambos puntos (verde y rojo) tienen los mismos valores de coordenadas en los ejes X e Y, en este caso las correspondientes a la ciudad de Toledo.
Sin embargo, presentan distintos CRS. Por este motivo, al representar ambos puntos en un mapa, se observa que no se están refiriendo a la misma localización geográfica. Esto es así porque el CRS define la referencia (punto x=0 e y =0) y las unidades de los ejes (grados, metros, millas).
Como conclusión, además de disponer de las coordenadas de los datos espaciales, es necesario conocer el CRS en el que están definidos para conocer de manera exacta su localización geográfica. Además, nótese que para cualquier análisis de datos espaciales es necesario que todos los geodatos se encuentren referenciados en el mismo CRS. Esto se consigue transformando (o proyectando) los datos a un CRS común, nunca sobreescribiendo el CRS de los mismos.
0.0.1 Tipos de CRS
A continuación se definen los dos grandes tipos de CRS, los CRS geográficos y los CRS proyectados.
0.0.1.1 CRS geográficos
Los CRS geográficos son aquellos en los que los parámetros empleados para localizar una posición espacial son la latitud y la longitud:
Latitud: Es la distancia angular expresada en grados sobre el plano definido por el ecuador terrestre. Determina la posición sobre de una localización en el eje Norte-Sur de la Tierra y toma valores en el rango \([-90º,90º]\) . Las líneas imaginarias determinadas por una sucesión de puntos con la misma latitud a lo largo del eje Este-Oeste se denominan paralelos. (Ver Fig. 2)
Longitud: Es la distancia angular expresada en grados sobre el plano definido por el meridiano de Greenwich. Determina la posición sobre de una localización en el eje Este-Oeste de la Tierra y toma valores en el rango \([-180º,180º]\) . Las líneas imaginarias determinadas por una sucesión de puntos con la misma longitud a lo largo del eje Este-Oeste se denominan meridianos (Ver Fig. 2)
Figure 2: Paralelos y Meridianos terrestres
Es muy importante destacar que en un sistema de coordenadas geográfico, es decir, basado en latitudes y longitudes, las distancias entre dos puntos representan distancias angulares. Por ejemplo, la distancia entre el meridiano de Greenwich y el meridiano correspondiente a la longitud 20º siempre es de +20º. Sin embargo, debido a la forma esférica de la Tierra, la longitud en metros entre ambos meridianos no es constante:
Figure 3: Distancia entre meridianos en distintas latitudes
0.0.1.2 CRS proyectados
La representación de formas tridimensionales en un soporte plano (dos dimensiones) presenta algunos retos. Por ello, es habitual trabajar con proyecciones de mapas.
Una proyección geográfica es un método para reducir la superficie de la esfera terrestre a un sistema cartesiano de dos dimensiones. Para ello, es necesario transformar las coordenadas longitud y latitud en coordenadas cartesianas x e y.
Es importante destacar que las proyecciones pueden incluir un punto de origen (X=0, Y=0) y unas unidades de distancia (habitualmente metros) específicas. Por ejemplo, la proyección cónica equiáreas de Albers (específica para Estados Unidos) define su punto de referencia (0,0) en la latitud 40º N y longitud 96º, y la unidad de variación están definida en metros. De ahí la importancia de conocer el CRS de los datos geográficos, como se expuso al principio de este tema.
Existen varias familias de proyecciones, que se pueden clasificar de diversas maneras: por tipo de superficie de proyección y por métrica a preservar.
a) Por tipo de superficie de proyección
El proceso de trasladar puntos de una esfera a un plano puede plantearse de manera práctica como el ejercicio de envolver una esfera con una superificie plana (como una hoja de papel) y trasladar los puntos de la esfera de manera lineal al punto de la superficie plana más cercano a ella. La Fig. 4 muestra estos tres tipos de proyección.
Figure 4: Tipos de proyección por superficie de proyección
A partir de este ejercicio, se plantean tres posibles soluciones, dependiendo del tipo de superficie que se use para proyectar.
- Proyecciones cilíndricas: Son aquellas proyecciones donde la superficie de proyección conforma un cilindro alrededor de la Tierra. Una de las proyecciones cilíndricas más conocidas es la proyección de Mercator (Ver Fig. 5).
Figure 5: Proyección Mercator
- Proyecciones cónicas: En este tipo de proyecciones, se plantea la superficie de proyección como una forma cónica. Como ejemplo, la proyección cónica equiáreas de Albers es una de las proyecciones que más suele usarse en la representación de mapas de América del Norte (Ver Fig. 6).
Figure 6: Proyección cónica equiáreas de Albers
- Proyecciones acimutales o planares: En este tipo de proyección se proyecta una porción de la Tierra sobre un plano que es tangente a la misma en el punto de referencia. Como ejemplos de proyecciones acimutales podemos destacar la proyección ortográfica (Ver Fig. 7).
Figure 7: Proyección ortogonal
b) Por métrica a preservar
Es importante tener en cuenta que cualquier proyección de la superficie de la Tierra produce distorsiones en una o varias características geográficas. Como ejemplos clásicos, la proyección de Mercator produce distorsiones del tamaño especialmente en aquellas regiones más cercanas a los polos (Groenlandia, que la proyección de Mercator presenta una área similar a la de África, tiene menor superficie real que Argelia). Otras de las métricas que suele verse distorsionada son la distancia entre dos puntos geográficos, la dirección o la forma de regiones de la Tierra.
A lo largo de la Historia se han desarrollado diversas proyecciones cuyo objetivo es preservar alguna o varias de las propiedades mencionadas anteriormente, sin embargo es importante destacar que no existe una proyección que sea capaz de preservar todas las métricas a la vez.
Según la metrica a presevar, las proyecciones se pueden clasificar en:
- Proyecciones conformales: Intentan preservar los ángulos que se forman en la superficie terrestre. Por ejemplo, la proyección de Mercator representa ángulos rectos en las intersecciones de los paralelos y los meridianos (Ver Fig. 8).
Figure 8: Ejemplo de proyección conformal: Mercator
- Proyecciones equivalentes: Preservan las proporciones de las áreas, provocando a su vez deformaciones en el resto de características, como la forma o los ángulos. La proyección acimutal equivalente de Lambers es un tipo de proyección equivalente (Ver Fig. 9).
Figure 9: Ejemplo de proyección equivalente: Proyección acimutal equivalente de Lambers
- Proyecciones equidistantes: Preservan la distancia entre dos puntos geográficos específicos. Por ejemplo, la proyección Plate carré preserva la distancia entre el Polo Norte y el Polo Sur (Ver Fig. 10).
Figure 10: Ejemplo de proyección equidistante: Platé carre
- Proyecciones de compromiso: No intentan preservar ninguna métrica en concreto. En su lugar, se centran en intentar encontrar un equilibrio entre las diversas distorsiones que provocan para intentar dar una representación más o menos representativa de la superficie terrestre. La proyección de Winkel Tripel, usada en los mapas de National Geographic, es un ejemplo de proyección de compromiso (Ver Fig. 11).
Figure 11: Ejemplo de proyección de compromiso: Winkel Tripel
En los anteriores ejemplos se ha añadido a cada proyección la indicatriz de Tissot. Ésta consiste en una serie de círculos imaginarios de igual área distribuidos sobre la superficie esférica de la Tierra en determinados puntos. De este manera, al presentar la indicatriz de Tissot en una proyección específica, se puede entender de una manera intuitiva la distorsión provocada por dicha proyección, ya que los círculos se ven distorsionados o preservados según los parámetros y la naturaleza de la proyección en cuestión.
0.0.2 Trabajando con proyecciones en R
Existe toda una serie de proyecciones predefinidas, identificadas mediante los códigos EPSG, ESRI, WKT o proj4 (en desuso en R, pero todavía admitidos). Existen varios recursos web donde se pueden consultar y seleccionar los códigos correspondientes:
Algunos de los códigos de proyecciones que es fundamental conocer son:
EPSG: 4326: Proyección correspondiente a WGS 84, que es el sistema usado por los sistemas GPS. Cuando trabajemos con coordenadas geográficas longitud/latitud, este es habitualmente el CRS de referencia.
EPSG: 3857: Código correspondiente a la proyección de Mercator, usada habitualmente por servicios como Google Maps, etc.
En la sección 0.0.3 veremos cómo encontrar un CRS usando el paquete
crsuggest.
El paquete sf permite obtener los parámetros de cualquier proyección mediante
la función st_crs():
library(sf)
# Ejemplo: EPSG WGS 84 (Sistema Global GPS): EPSG 4326
st_crs(4326)
#> Coordinate Reference System:
#> User input: EPSG:4326
#> wkt:
#> GEOGCRS["WGS 84",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> CS[ellipsoidal,2],
#> AXIS["geodetic latitude (Lat)",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["geodetic longitude (Lon)",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> USAGE[
#> SCOPE["Horizontal component of 3D system."],
#> AREA["World."],
#> BBOX[-90,-180,90,180]],
#> ID["EPSG",4326]]
# Usando código ESRI North America Albers Equal Area Conic
st_crs("ESRI:102008")
#> Coordinate Reference System:
#> User input: ESRI:102008
#> wkt:
#> PROJCRS["North_America_Albers_Equal_Area_Conic",
#> BASEGEOGCRS["NAD83",
#> DATUM["North American Datum 1983",
#> ELLIPSOID["GRS 1980",6378137,298.257222101,
#> LENGTHUNIT["metre",1]]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["Degree",0.0174532925199433]]],
#> CONVERSION["North_America_Albers_Equal_Area_Conic",
#> METHOD["Albers Equal Area",
#> ID["EPSG",9822]],
#> PARAMETER["Latitude of false origin",40,
#> ANGLEUNIT["Degree",0.0174532925199433],
#> ID["EPSG",8821]],
#> PARAMETER["Longitude of false origin",-96,
#> ANGLEUNIT["Degree",0.0174532925199433],
#> ID["EPSG",8822]],
#> PARAMETER["Latitude of 1st standard parallel",20,
#> ANGLEUNIT["Degree",0.0174532925199433],
#> ID["EPSG",8823]],
#> PARAMETER["Latitude of 2nd standard parallel",60,
#> ANGLEUNIT["Degree",0.0174532925199433],
#> ID["EPSG",8824]],
#> PARAMETER["Easting at false origin",0,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8826]],
#> PARAMETER["Northing at false origin",0,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8827]]],
#> CS[Cartesian,2],
#> AXIS["(E)",east,
#> ORDER[1],
#> LENGTHUNIT["metre",1]],
#> AXIS["(N)",north,
#> ORDER[2],
#> LENGTHUNIT["metre",1]],
#> USAGE[
#> SCOPE["Not known."],
#> AREA["North America - onshore and offshore: Canada - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon. United States (USA) - Alabama; Alaska (mainland); Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming."],
#> BBOX[23.81,-172.54,86.46,-47.74]],
#> ID["ESRI",102008]]
# Usando proj4string: Robinson
st_crs("+proj=robin")
#> Coordinate Reference System:
#> User input: +proj=robin
#> wkt:
#> PROJCRS["unknown",
#> BASEGEOGCRS["unknown",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ID["EPSG",6326]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8901]]],
#> CONVERSION["unknown",
#> METHOD["Robinson"],
#> PARAMETER["Longitude of natural origin",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8802]],
#> PARAMETER["False easting",0,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8806]],
#> PARAMETER["False northing",0,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8807]]],
#> CS[Cartesian,2],
#> AXIS["(E)",east,
#> ORDER[1],
#> LENGTHUNIT["metre",1,
#> ID["EPSG",9001]]],
#> AXIS["(N)",north,
#> ORDER[2],
#> LENGTHUNIT["metre",1,
#> ID["EPSG",9001]]]]La mayoría de los objetos espaciales serán de la clase sf, por tanto, resulta
interesante conocer cómo se proyectan estos objetos.
Es posible proyectar un objeto sf mediante la función st_transform(). En el
siguiente ejemplo vemos cómo partimos de un objeto con EPSG:4326 y cambiamos
su proyección a otras proyecciones, como Mercator o Robinson (Fig.
??).
# Usa datos del paquete mapSpain
library(giscoR)
paises <- gisco_get_countries()
# Comprobamos el CRS de estos datos
# Se puede almacenar en un objeto y usar posteriormente
# Vemos que es EPSG:4326, por tanto son coordenadas geográficas longitud/latitud
st_crs(paises)
#> Coordinate Reference System:
#> User input: EPSG:4326
#> wkt:
#> GEOGCS["WGS 84",
#> DATUM["WGS_1984",
#> SPHEROID["WGS 84",6378137,298.257223563,
#> AUTHORITY["EPSG","7030"]],
#> AUTHORITY["EPSG","6326"]],
#> PRIMEM["Greenwich",0,
#> AUTHORITY["EPSG","8901"]],
#> UNIT["degree",0.0174532925199433,
#> AUTHORITY["EPSG","9122"]],
#> AUTHORITY["EPSG","4326"]]
# Plot
plot(st_geometry(paises), axes = TRUE)
Figure 12: Proyección del mundo en coordenadas geográficas (EPSG 4326)
# Proyectamos a Mercator
# El eje cambia porque Mercator usa metros
paises_merc <- st_transform(paises, st_crs(3857))
plot(st_geometry(paises_merc), axes = TRUE)
Figure 13: Proyección del mundo en Mercator (EPSG 3857)
# Proyectamos a Robinson
paises_robin <- st_transform(paises, st_crs("+proj=robin"))
plot(st_geometry(paises_robin), axes = TRUE)
Figure 14: Proyección del mundo en Robinson (+proj=robin)
Como se comentó anteriormente, cuando se usan geodatos de diversas fuentes, es necesario que todos presenten el mismo CRS. En la Fig ?? se muestra lo que ocurre si esto no se cumple:
# Añadimos a este mapa puertos mundiales de giscoR
puertos <- gisco_get_ports()
plot(st_geometry(paises_robin), main = "Puertos en el mundo")
plot(st_geometry(puertos), add = TRUE, col = "red", pch = 20)
Figure 15: Ejemplo: Puertos del mundo, CRS no alineados
# Ha habido algun error... Comprueba CRS
st_crs(puertos) == st_crs(paises_robin)
#> [1] FALSE
# Los puertos no están en Robinson! Proyectamos al mismo CRS
puertos_robin <- st_transform(puertos, st_crs(paises_robin))
plot(st_geometry(paises_robin), main = "Puertos en el mundo")
plot(st_geometry(puertos_robin), add = TRUE, col = "blue", pch = 20)
Figure 16: Ejemplo: Puertos del mundo, CRS alineados
Como vemos, en el primer mapa los puertos se concentran en un único punto, dado que no están referenciados en el mismo CRS. Tras proyectarlos, el mapa se representa adecuadamente.
En otros paquetes, como sp o raster, existen funciones parecidas. Cuando
empleemos el paquete sp podemos usar las funciones CRS() y spTransform():
library(sp)
# Convertimos sf a sp
paises_sp <- as(paises, "Spatial")
# En sp podemos usar:
# CRS("+proj=robin")
#
# O también desde sf
# CRS(st_crs(paises_robin)$proj4string)
paises_sp_robin <- spTransform(paises_sp, CRS("+proj=robin"))
plot(paises_sp_robin)
Figure 17: Transformaciones en sp
En el caso de un objeto raster, podemos usar crs() y projectRaster():
library(raster)
# Extrae información de altitud para Paises Bajos
elev <- getData("alt", country = "NLD", path = tempdir())
# Transforma
elev_robinson <- projectRaster(elev, crs = crs("+proj=robin"))
plot(elev_robinson)
Figure 18: Transformaciones en raster
Por último, en el paquete terra las funciones correspondientes son crs() y
project():
library(terra)
# Convierte de raster a terra
elev_terra <- rast(elev)
# Transforma
elev_terra_robinson <- terra::project(elev_terra, terra::crs(elev_terra))
plot(elev_terra_robinson)
Figure 19: Transformaciones en terra
0.0.3 ¿Qué proyección uso?
El CRS adecuado para cada análisis depende de la localización y el rango
espacial de los datos. Un CRS adecuado para representar un mapa del mundo puede
no serlo para representar datos de zonas específicas de la Tierra. Los recursos
web mencionados anteriormente permiten la búsqueda de CRS por zona geográfica, y
adicionalmente en R existe el paquete crsuggest (Walker, 2021) que nos
facilita la labor, sugiriendo el CRS más adecuado para cada zona:
library(crsuggest)
# CRS para Paises Bajos
# Usando raster
sugerencias <- suggest_crs(elev)| crs_code | crs_name | crs_type | crs_gcs | crs_units | crs_proj4 |
|---|---|---|---|---|---|
| 28992 | Amersfoort / RD New | projected | 4289 | m | +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812 +units=m +no_defs |
| 28991 | Amersfoort / RD Old | projected | 4289 | m | +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=0 +y_0=0 +ellps=bessel +towgs84=565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812 +units=m +no_defs |
| 5651 | ETRS89 / UTM zone 31N (N-zE) | projected | 4258 | m | +proj=tmerc +lat_0=0 +lon_0=3 +k=0.9996 +x_0=31500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs |
| 5649 | ETRS89 / UTM zone 31N (zE-N) | projected | 4258 | m | +proj=tmerc +lat_0=0 +lon_0=3 +k=0.9996 +x_0=31500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs |
| 3812 | ETRS89 / Belgian Lambert 2008 | projected | 4258 | m | +proj=lcc +lat_0=50.797815 +lon_0=4.35921583333333 +lat_1=49.8333333333333 +lat_2=51.1666666666667 +x_0=649328 +y_0=665262 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs |
| 3447 | ETRS89 / Belgian Lambert 2005 | projected | 4258 | m | +proj=lcc +lat_0=50.797815 +lon_0=4.35921583333333 +lat_1=49.8333333333333 +lat_2=51.1666666666667 +x_0=150328 +y_0=166262 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs |
| 31370 | Belge 1972 / Belgian Lambert 72 | projected | 4313 | m | +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-99.059,53.322,-112.486,0.419,-0.83,1.885,-1 +units=m +no_defs |
| 31300 | Belge 1972 / Belge Lambert 72 | projected | 4313 | m | +proj=lcc +lat_0=90 +lon_0=4.35693972222222 +lat_1=49.8333333333333 +lat_2=51.1666666666667 +x_0=150000.01256 +y_0=5400088.4378 +ellps=intl +towgs84=-99.059,53.322,-112.486,0.419,-0.83,1.885,-1 +units=m +no_defs |
| 21500 | Belge 1950 (Brussels) / Belge Lambert 50 | projected | 4809 | m | +proj=lcc +lat_0=90 +lon_0=0 +lat_1=49.8333333333333 +lat_2=51.1666666666667 +x_0=150000 +y_0=5400000 +ellps=intl +pm=brussels +units=m +no_defs |
| 23095 | ED50 / TM 5 NE | projected | 4230 | m | +proj=tmerc +lat_0=0 +lon_0=5 +k=0.9996 +x_0=500000 +y_0=0 +ellps=intl +towgs84=-89.5,-93.8,-123.1,0,0,-0.156,1.2 +units=m +no_defs |
# Probamos sugerencia
crs_suggest <- suggest_crs(elev, limit = 1)
elev_suggest <- projectRaster(elev, crs = raster::crs(crs_suggest$crs_proj4))
plot(elev_suggest)
Figure 20: raster: Ejemplo de transformación usando crsuggest
# Ejemplo con sf: China
china <- gisco_get_countries(country = "China")
china_crs <- suggest_crs(china, limit = 1)
china_suggest <- st_transform(
china,
st_crs(as.integer(china_crs$crs_code))
)
plot(st_geometry(china_suggest), axes = TRUE)
Figure 21: sf: Ejemplo de transformación usando crsuggest